packages<-c('ggplot2','dplyr','mada')
lapply(packages, library, character.only=T)
## [[1]]
## [1] "ggplot2" "stats" "graphics" "grDevices" "utils" "datasets"
## [7] "methods" "base"
##
## [[2]]
## [1] "dplyr" "ggplot2" "stats" "graphics" "grDevices" "utils"
## [7] "datasets" "methods" "base"
##
## [[3]]
## [1] "mada" "mvmeta" "ellipse" "mvtnorm" "dplyr"
## [6] "ggplot2" "stats" "graphics" "grDevices" "utils"
## [11] "datasets" "methods" "base"
"%&%" = function(a,b) paste(a,b,sep="")
h2.dir <- '/Volumes/im-lab/nas40t2/Data/Annotations/heritability/'
gcta <- read.table('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/expArch_DGN-WB_imputedGTs/DGN-WB.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-13.txt',header=TRUE) %>% dplyr::select(gene,local.h2)
dim(gcta)
## [1] 12719 2
##get price h2 (custom-made human array containing 23,720 unique oligonucleotide probes, http://www.nature.com/nature/journal/v452/n7186/full/nature06758.html)
price <- read.table(h2.dir %&% "Alkes/h2all.txt",header=TRUE) %>% dplyr::rename(gene=gname)
dim(price)
## [1] 15078 7
wright <- read.table(h2.dir %&% "Fred/h2PB.txt",header=TRUE) %>% dplyr::rename(gene=gname)
dim(wright)
## [1] 17792 2
all<-full_join(gcta,price,by='gene')
all<-full_join(all,wright,by='gene')
#add GTEx whole blood
gtex <- read.table('/Volumes/im-lab/nas40t2/hwheeler/cross-tissue/GTEx.TW.WholeBlood.h2.all.models_FHSfdr0.05.Chr1-22_globaleQTLOtherChr_reml-no-constrain.2015-12-14.txt',header=T) %>% mutate(gtexWB=local.h2) %>% dplyr::select(gene,gtexWB)
all<-full_join(all,gtex,by='gene')
dim(all)
## [1] 22567 10
nona <- all[complete.cases(all),]
dim(nona)
## [1] 5518 10
#dgn v. gtex whole blood
ggplot(nona,aes(x=local.h2,y=gtexWB))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()+coord_cartesian(ylim=c(-0.15,1.02))
res<-cor.test(nona$local.h2,nona$gtexWB,method='s')
CIrho(res$estimate,length(nona$local.h2))
## rho 2.5 % 97.5 %
## [1,] 0.2713322 0.2467124 0.2956019
#dgn v. price blood
ggplot(nona,aes(x=local.h2,y=h2bloodcis))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()
res<-cor.test(nona$local.h2,nona$h2bloodcis,method='s')
CIrho(res$estimate,length(nona$local.h2))
## rho 2.5 % 97.5 %
## [1,] 0.3245192 0.3007081 0.3479261
#dgn v. wright peripheral blood
ggplot(nona,aes(x=local.h2,y=h2PB))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()
res<-cor.test(nona$local.h2,nona$h2PB,method='s')
CIrho(res$estimate,length(nona$local.h2))
## rho 2.5 % 97.5 %
## [1,] 0.07345576 0.04716111 0.09964868
#wright v. price blood
ggplot(nona,aes(x=h2PB,y=h2bloodcis))+geom_point(alpha=1/4)+geom_smooth()+theme_bw()
res<-cor.test(nona$h2PB,nona$h2bloodcis,method='s')
CIrho(res$estimate,length(nona$local.h2))
## rho 2.5 % 97.5 %
## [1,] 0.1064164 0.08025564 0.1324306
plot all 3 blood
nona <- all[complete.cases(all),]
dim(nona)
## [1] 5518 10
#dgn v. price blood, color by wright
ggplot(nona,aes(x=local.h2,y=h2bloodcis,col=h2PB))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")
#dgn v. wright peripheral blood, color by price
ggplot(nona,aes(x=local.h2,y=h2PB,col=h2bloodcis))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")
#wright v. price blood, color by dgn
ggplot(nona,aes(x=h2PB,y=h2bloodcis,col=local.h2))+geom_point()+geom_smooth()+theme_bw()+ scale_colour_gradient(low="pink",high="blue")
cor tables
round(cor(all[,-1],use='pairwise.complete.obs',method="s"),3)
## local.h2 h2adip h2blood h2adipcis h2adiptra h2bloodcis
## local.h2 1.000 0.126 0.317 0.134 -0.006 0.307
## h2adip 0.126 1.000 0.335 0.328 0.617 0.231
## h2blood 0.317 0.335 1.000 0.246 0.140 0.392
## h2adipcis 0.134 0.328 0.246 1.000 -0.415 0.235
## h2adiptra -0.006 0.617 0.140 -0.415 1.000 0.059
## h2bloodcis 0.307 0.231 0.392 0.235 0.059 1.000
## h2bloodtra 0.041 0.172 0.629 0.089 0.138 -0.339
## h2PB 0.083 0.111 0.190 0.066 0.074 0.124
## gtexWB 0.281 0.051 0.110 0.060 -0.007 0.120
## h2bloodtra h2PB gtexWB
## local.h2 0.041 0.083 0.281
## h2adip 0.172 0.111 0.051
## h2blood 0.629 0.190 0.110
## h2adipcis 0.089 0.066 0.060
## h2adiptra 0.138 0.074 -0.007
## h2bloodcis -0.339 0.124 0.120
## h2bloodtra 1.000 0.102 0.004
## h2PB 0.102 1.000 0.042
## gtexWB 0.004 0.042 1.000
#Blood only
allblood <- dplyr::mutate(all,dgnWB=local.h2,Price=h2bloodcis,Wright=h2PB) %>% dplyr::select(dgnWB,gtexWB,Price,Wright)
#Pearson
round(cor(allblood,use='pairwise.complete.obs'),3)
## dgnWB gtexWB Price Wright
## dgnWB 1.000 0.469 0.391 0.150
## gtexWB 0.469 1.000 0.249 0.118
## Price 0.391 0.249 1.000 0.156
## Wright 0.150 0.118 0.156 1.000
#Spearman
round(cor(allblood,use='pairwise.complete.obs',method="s"),3)
## dgnWB gtexWB Price Wright
## dgnWB 1.000 0.281 0.307 0.083
## gtexWB 0.281 1.000 0.120 0.042
## Price 0.307 0.120 1.000 0.124
## Wright 0.083 0.042 0.124 1.000
cor tables - intersection of the 4 whole blood studies
dim(nona)
## [1] 5518 10
round(cor(nona[,-1],use='pairwise.complete.obs',method="s"),3)
## local.h2 h2adip h2blood h2adipcis h2adiptra h2bloodcis
## local.h2 1.000 0.113 0.326 0.125 -0.012 0.325
## h2adip 0.113 1.000 0.272 0.374 0.572 0.234
## h2blood 0.326 0.272 1.000 0.261 0.065 0.450
## h2adipcis 0.125 0.374 0.261 1.000 -0.409 0.275
## h2adiptra -0.012 0.572 0.065 -0.409 1.000 0.019
## h2bloodcis 0.325 0.234 0.450 0.275 0.019 1.000
## h2bloodtra 0.026 0.107 0.580 0.064 0.105 -0.337
## h2PB 0.073 0.050 0.147 0.041 0.041 0.106
## gtexWB 0.271 0.041 0.087 0.057 -0.023 0.130
## h2bloodtra h2PB gtexWB
## local.h2 0.026 0.073 0.271
## h2adip 0.107 0.050 0.041
## h2blood 0.580 0.147 0.087
## h2adipcis 0.064 0.041 0.057
## h2adiptra 0.105 0.041 -0.023
## h2bloodcis -0.337 0.106 0.130
## h2bloodtra 1.000 0.066 -0.035
## h2PB 0.066 1.000 0.032
## gtexWB -0.035 0.032 1.000
#Blood only
allblood <- dplyr::mutate(nona,dgnWB=local.h2,Price=h2bloodcis,Wright=h2PB) %>% dplyr::select(dgnWB,gtexWB,Price,Wright)
#Pearson
round(cor(allblood,use='pairwise.complete.obs'),3)
## dgnWB gtexWB Price Wright
## dgnWB 1.000 0.461 0.407 0.134
## gtexWB 0.461 1.000 0.255 0.111
## Price 0.407 0.255 1.000 0.149
## Wright 0.134 0.111 0.149 1.000
#Spearman
round(cor(allblood,use='pairwise.complete.obs',method="s"),3)
## dgnWB gtexWB Price Wright
## dgnWB 1.000 0.271 0.325 0.073
## gtexWB 0.271 1.000 0.130 0.032
## Price 0.325 0.130 1.000 0.106
## Wright 0.073 0.032 0.106 1.000